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Wilson's numerical renormalization group (NRG) method 
for the calculation of dynamic properties of impurity models is 
generalized to investigate the effective impurity model of the 
dynamical mean field theory at finite temperatures. We cal- 
culate the spectral function and self-energy for the Hubbard 
model on a Bethe lattice with infinite coordination number di- 
rectly on the real frequency axis and investigate the phase dia- 
gram for the Mott-Hubbard metal-insulator transition. While 
for T < T c ~ 0.02W (W: bandwidth) we find hysteresis with 
first-order transitions both at U c i (defining the insulator to 
metal transition) and at U C 2 (defining the metal to insula- 
tor transition), at T > T c there is a smooth crossover from 
metallic-like to insulating-like solutions. 

7f.10.Fd, 71.30.+h 



I. INTRODUCTION 

During the past decade, the development and appli- 
cation of the dynamical mean field theory (DMFT) has 
led to a considerable increase in our understanding of 
strongly correlated electron systems. The DMFT has 
originally been derived from the limit of infinite spatial 
dimensionality (or, equivalently, infinite lattice connec- 
tivity) of lattice fermion models, such as the Hubbard 
model Jjl[. In this limit, the self-energy becomes purely 
local [g]7 which is a consequence of the required scaling 
of the hopping matrix element t — t* j 1 \fd, with t* fixed 
and d the lattice dimension. 

It has been realized in the work of Jarrell || and 
Georges and Kotliar Q that such a local self-energy 
can be calculated from a much simpler, but nevertheless 
highly non-trivial model: the single impurity Anderson 
model (SIAM) @. The self-energy of the SIAM is local 
because the Coulomb correlation in this model only acts 
on the impurity site. The difference between the SIAM 
and the lattice model under consideration is then built 
in via a self-consistency condition [0]. In this way, the 
DMFT became a powerful tool for the investigation of 
various lattice models such as the Hubbard model and 
the periodic Anderson model (for a review see ||). The 
success of this approach, however, depends on the avail- 
ability of reliable methods for the calculation of the self- 
energy of an effective SIAM. Perturbative methods, such 
as the iterated perturbation theory || or the non-crossing 
approximation H , have been shown to give qualitatively 



correct results for a variety of physical problems. The 
numerical implementation of these methods allows one 
to solve the impurity model with a minimum of compu- 
tational effort (typically a few seconds on a workstation) 
so that the relevant parameter space of the model can be 
scanned very quickly. 

However, most of the phenomena of interest in strongly 
correlated systems are inherently non-perturbative, so 
that none of the parameters in the Hamiltonian can be re- 
garded as a small perturbation. In general we therefore 
have to apply non-perturbative methods, even in cases 
where perturbative approaches such as the iterated per- 
turbation theory or the non-crossing approximation ap- 
pear to give a complete picture of the solution. 

The most widely used non-perturbative method in this 
context is the quantum Monte-Carlo approach |||| . The 
advantages of this method are its flexibility (a wide range 
of physical problems can be studied with only relatively 
minor changes in the program) and the possibility of ob- 
taining a 'numerically exact' solution of G(r), the single- 
particle Green function on the imaginary time axis. The 
main disadvantage of the quantum Monte-Carlo method 
is the drastic increase of computation time upon either 
increasing the Coulomb repulsion U or decreasing the 
temperature T. Furthermore, the analytic continuation 
of the data on the imaginary time or frequency axis to 
the real axis represents a difficult and numerically ill- 
conditioned problem (see || for the application of the 
maximum entropy method to this problem). 

Another non-perturbative method applicable here is 
the exact diagonalization technique (see, e.g., P,p|,^0t). 
In this method, the continuous conduction band of the 
effective SIAM is approximated by a discrete set of states 
(approximately 8-12 states). The value of U does not 
impose a problem here as the impurity (together with the 
conduction electron states) is diagonalized exactly. The 
main disadvantage of the exact diagonalization technique 
is its inability to resolve low energy features such as a 
narrow quasiparticle resonance at the Fermi level. 

The above-mentioned restrictions concerning the value 
of U and T, or the low energy resolution, do not apply 
to the numerical renormalization group method (NRG) 
PUP^I which has only recently been used to investigate 
lattice models within the DMFT JTfHT7|- The NRG as 
well has its drawbacks, which will be discussed in Sec. II 
of this paper; nevertheless one would expect the NRG 
to be an ideal tool to calculate the self-energy of the ef- 
fective Anderson model in the DMFT, simply because 
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it has proven to be very successful in the investigation 
of the physics of the standard SIAM. For example, the 
NRG method is able to resolve, both in static and dy- 
namic properties, the exponentially small Kondo scale 
for large values of U (which can be seen neither in quan- 
tum Monte-Carlo nor in exact diagonalization). One can 
also study in detail the scaling spectrum of the quasi- 
particle peak (T^Jl^] and the temperature dependence of 
transport properties @|IJ of the SIAM. 

Applications of the NRG within the DMFT include the 
investigation of the Mott-transition |H| , the problem 
of charge ordering in the extended Hubbard model [jl6| , 
and the formation of the heavy-fermion liquid in the pe- 
riodic Anderson model |17j] . In all these investigations, 
the temperature was restricted to T = 0. 

In this paper, we present a study of a strongly cor- 
related lattice model within the DMFT by applying the 
NRG method at finite temperatures [ pOpl] ]. In particu- 
lar, we address a problem which has been the topic of an 
intense debate over the last couple of years: the details 
of the Mott-transition from a paramagnetic metal to a 
paramagnetic insulator in the half-filled Hubbard model 

|[ljlH§. 

The paper is organized as follows: the NRG method 
is introduced in Sec. II, with particular emphasis on the 
calculation of finite temperature dynamics. In Sec. Ill, 
the previous results for the Mott-transition in the Hub- 
bard model (within DMFT) are discussed. The results 
from the NRG for the finite temperature Mott-transition 
are then presented in Sec. IV. The paper is concluded 
with a summary in Sec. V. 

II. THE NUMERICAL RENORMALIZATION 
GROUP METHOD AT FINITE TEMPERATURES 

A. General Concepts 

The basic ideas of the NRG method were developed 
by Wilson for the investigation of the Kondo model . 
Krishna- murthy, Wilkins, and Wilson |l2| later applied 
the NRG to a related model, the single impurity Ander- 
son model (SIAM) with the Hamiltonian 

(J 

+E^ c L^ + E^(/^ + c La) • (i) 

ko ka 

In the model (Q), cj^ denote annihilation (creation) op- 
erators for band states with spin a and energy f„ 
those for impurity states with spin a and energy S{. The 
Coulomb interaction for two electrons at the impurity 
site is given by U and both subsystems are coupled via 
a hybridization V. 

The hybridization function 



k 

is usually assumed to be constant between the band-edges 
—D and D, but will acquire some frequency dependence 
in the effective Anderson model within the DMFT (the 
necessary changes in the NRG-procedure due to the non- 
constant A(lu) were discussed in 

The first step to set up the rcnormalization group 
transformation is a logarithmic discretization of the con- 
duction band: the continuous conduction band is divided 
into infinitely many intervals [£ n +i>£n] an d [— — £n+i] 
with £ n = DA~ n and n = 0, 1, 2, . . . , oo. Here, A is 
the NRG discretization parameter (typical values used 
in the calculations are A = 1.5, . . . , 2). The conduction 
band states in each interval are then replaced by a sin- 
gle state. While this approximation by a discrete set of 
states involves some coarse graining at higher energies, it 
captures arbitrarily small energies near the Fermi level. 

In a second step, this discrete model is mapped on a 
semi-infinite chain form via a tridiagonalization proce- 
dure (for details, see @Q and section 4.2 in Q). The 
Hamiltonian of the semi-infinite chain has the following 
form: 

H = J2 *fti* + Uf]hf\h + E v(f$<%° + cL/ CT ) 

oo 

+ tn fc^o-Cn+lCT + c n+-LtT C n<rJ ■ (3) 

o",n— 

This form is valid for a general symmetric conduction 
band density of states. The impurity now couples only 
to a single fermionic degree of freedom (the Cgj^), with a 
hybridization V. Due to the logarithmic discretization, 
the hopping matrix elements decrease as t n cx A~"/ 2 . 
This means that, in going along the chain, the parameters 
in the Hamiltonian evolve from high energies (given by 
D and U) to arbitrarily low energies (given by DA""/ 2 ). 
The renormalization group transformation is now set up 
in the following way. 

We start with the solution of the isolated impurity, 
i.e., with the knowledge of all eigenstates, eigenenergies, 
and matrix elements. The first step of the renormal- 
ization group transformation is to add the first conduc- 
tion electron site, set up the Hamiltonian matrices for 
the enlarged Hilbert space, and obtain the new eigen- 
states, eigenenergies, and matrix elements by diagonal- 
izing these matrices. This procedure is then iterated. 
An obvious problem occurs after only a few steps of the 
iteration. The Hilbert space grows as 4 N (with N the 
size of the cluster) , which makes it impossible to keep all 
the states in the calculation. Wilson therefore devised 
a very simple truncation procedure in which only those 
states with the lowest energies (typically a few hundred) 
are kept. This truncation scheme is very successful but 
relies on the fact that the hopping matrix elements are 
falling off exponentially. High-energy states therefore do 
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not change the very low frequency behaviour and can be 
neglected. This procedure gives for each cluster a set of 
eigenenergies and matrix elements from which a number 
of physical properties can be derived. 



B. Finite Temperature Dynamics 

Here we want to discuss in detail the calculation of the 
finite temperature spectral function 



with 



AJlo) = --Im G a (io + iS + ) 

7T 



G a (z) = i dt e^{[f a (t)Jt] + ) . 



(4) 
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From the iterative diagonalization described above, one 
can easily calculate the spectral functions for each cluster 
of size N via lEl 
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Here {j^^} and {jwi)^} are sets of eigenstates of the 
Hamiltonian for the cluster of size TV, and E^ are the 
corresponding eigenenergies and Zn the grand canonical 
partition function (the spin index a will be dropped in 
the following). As the length of the cluster is succes- 
sively increased, and the A a ^(uj) calculated in each step, 
Eq. (|^) defines a whole set of spectral functions. These 
data are combined to give spectral functions as shown, 
e.g., in Fig. || in the following way. 

Let us first describe theprocedure for calculating the 
T = spectral function [^,^0|. The diagonalization of 
the clusters N = 0,1,2,... yields the excitation spec- 
trum u) nm = E^ — E^ on a set of decreasing energy 
scales wo > co\ > ^2 > • ■ ■ (wjv is the smallest scale in 
the truncated Hamiltonian Hn, i.e., lon — tN and for 
a flat band one has lon ~ DA - ^ -1 )' 2 ). Excitations 
w <C lon are not described within cluster N. They are 
obtained accurately in subsequent iterations from larger 
clusters. Similarly, excitations u> u>n are outside the 
energy window for cluster N (whose width is limited on 
the high energy side by the truncation of the spectrum) . 
Information on these excitations is contained in previous 
iterations for some smaller cluster N' < N. It is therefore 
possible to use Eq. (||) for each N = 0, 1, . . . to calculate 
the T = spectral density at an appropriate set of de- 
creasing frequencies for each cluster. These frequencies 
are chosen to be lo ~ 2lon within the energy window of 
the cluster under consideration (whose width, in units of 
ujn, typically lies in the range 6-10 for A = 1.5 — 2.0). 

At finite temperature, the above procedure is modified 
as follows. For a given temperature T, which we identify 



with Tm ~ &m f° r some M , one evaluates the spectral 
density in Eq. (||) at the same characteristic frequencies 
u> = 2u>n as those used for the T = calculation, down 
to a minimum frequency corresponding to lo as T = Tm- 
Compared to the T — calculation, many more exci- 
tations will contribute at finite T, as shown in Fig. 0. 
When lo = Ilon becomes comparable to or smaller than 
the temperature of interest, T = Tm, it is clear that exci- 
tations will start to contribute to the spectral density at 
frequency lo which are not contained in cluster N . It is 
still possible to calculate the spectral density at frequen- 
cies lo — 2lon such that —Tm < w < +Tm by using the 
cluster of size M corresponding to the temperature. This 
is achieved by broadening the ^-functions in the spec- 
trum of cluster M with broadening functions of width T 
(see below). This gives a very good estimate of the lead- 
ing contribution to the spectral density for all frequencies 
M < T. It recovers, for example, the known Fermi liquid 
relations for transport quantities of the Anderson model 
p0| , pl| . Due to the limited resolution, proportional to 
Tm, the above scheme will, however, tend to broaden the 
spectral densities too much at higher temperatures. 
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FIG. 1. The spectrum of many-body excitations measured 
with respect to the ground-state energy E g — 0, and the 
possible transitions contributing to the single-particle spec- 
tral function, a) For T = 0, only transitions between the 
ground-state and excited states are possible; b) for T > 0, 
transitions between excited states are possible as well. The 
dotted line indicates the cut-off in the spectrum due to the 
truncation of states. 

This is not a problem for the finite temperature spec- 
tral densities presented in this paper. The reason, as we 
shall see below, is that the width of the Kondo resonance 
in the effective impurity model is always very much larger 
than the temperatures of interest (typically by a factor 
of 10 larger). 

The above scheme becomes increasingly accurate as 
the temperature is lowered, eventually connecting contin- 
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uously the finite and zero temperature spectral densities 
as T — > 0. 

There are several ways to put together the discrete in- 
formation from the clusters in order to arrive at contin- 
uous curves for spectral densities. One approach [ p0|j2l| 
replaces the delta functions in Eq. by appropriate 
broadening functions (see Eq. (Q-|8|)) and evaluates the 
spectral densities at the characteristic frequencies defined 
above. It is also possible to first combine information 
on the discrete spectra from successive clusters (N and 
N + 2, to avoid even/odd effects) and then broaden the 
spectra Q. Below, we describe this latter approach, 
which we used to obtain most results in this paper. A 
comparison between the two approaches gave only minor 
differences in the results for the spectral function. 
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FIG. 2. Superposition of the <5-peaks in the spectral density 
of all clusters up to length N [see (a)] with the 5-peaks of the 
cluster of length 7V+2 [see (b)] . This procedure gives the spec- 
tral information contained in all clusters up to length TV + 2 
[see (c)]. The spikes indicate the weight of the 8- functions in 
the spectral density, and the lines in (a) and (b) correspond to 
the weighting function as described in the text. The <5-peaks 
in the interval [w^ 2 ,Wmin] and for u) > w^i 2 appearing in 
(c) are identical to those appearing in (a) and (b), repectively, 
as indicated by the arrows. 

The starting point is the set of (5-peaks obtained for a 
small cluster of size N where the truncation is not yet 
effective (see Fig. ||a). The spectral distribution for the 
cluster of length (N + 2) is shown in Fig. ||b. The mini- 
mal frequency appearing in the spectrum for the cluster 



of length (N + 2), w mi „ , is reduced approximately by 
a factor of A compared to the frequency uJ^ in , while the 
maximal frequency io^ 2 is now determined by the num- 
ber of states retained after truncation. From the two sets 
of 5-peaks we keep those peaks which are in the interval 
[ w min 2 > w min] an d above uj^ 2 . The peaks in the overlap- 
ping region [w^ in ,Wmax 2 ] are taken from both the previ- 



ous clusters and the one of length N + 2, and are added 
with a weighting function which is, for simplicity, just 
a linear function with values from to 1 for arguments 
between u>^ in and oj^ 2 (for the previous clusters) and 
with values from 1 to (for the cluster of length N + 2) 
The resulting set of <5-peaks is shown in Fig. and 
can then be used to further iterate this procedure (with 
the cluster of length TV + 4, and so on), up to the cluster 
of length M defined by T = T M . 

The resulting spectrum is still discrete. To visualize 
the distribution of spectral weight it is convenient to 
broaden the J-peaks using appropriate broadening func- 
tions. For the results shown in this paper we used a 
Lorentzian 



S((jJ- LO n ) 



1 



2tt (w - uj n ) 2 + b 2 



(7) 



with width b = 0.6T for u) r , 
logarithmic scale 



< 4T and a Gaussian on a 



S(oj 



-b 2 /4 



■. exp 



(In u> — In uj n y 
P 



(8) 



with width b = 0.3 for u n > 4T. 

So far, we have not made any reference to the appli- 
cation of the NRG to the effective Anderson model in 
the DMFT. The necessary steps are described in [Q for 
the case of T = and can be used for finite tempera- 
tures equally well. In particular, the expression of the 
self-energy via 



£«,(<«;) = U 



(9) 



with the correlation function F a (u>) = ((.faftfa,fl))ui 
holds for both T — and T > (for a discussion of the 
advantage of using Eq. @ for the calculation of S(w), 
see@). 

Let us now comment on the choice of temperatures 
used in the calculations reported in this paper. It is clear 
from the above discussion that the temperatures are cho- 
sen to lie within the excitation spectrum of the cluster M 
for which the NRG iteration is terminated. Keeping the 
position of Tm within the excitation spectrum constant, 
one has to reduce Tm by a factor A when the largest 
cluster is of length M + 2. This defines a discrete set of 
temperatures Tjy = Tm^ n ^ m ^ 2 for which we perform 
the NRG calculations. 

For a variety of applications within the DMFT one 
would certainly prefer to vary the temperature continu- 
ously (to find, e.g., the critical temperatures for a phase 
transition). Such a continuous variation is difficult within 
NRG. It is certainly possible to achieve a large variation 
in temperature by a modest variation in A and using a 
fixed length of the cluster (due to the exponential de- 
pendence of Tm on A) . The results obtained in this way 
would, however, contain different systematic errors, as 
the accuracy of the NRG is enhanced upon reducing A. 
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One should therefore try to correct this A-dependence, 
e.g., by extrapolating the results to A — 1. We have not 
attempted to correct for the A-dependence and instead 
worked with a fixed A = 1.64 and different cluster sizes 
(the number of states retained after truncation is 600, 
not counting degeneracies). 



III. THE MOTT-HUBBARD METAL-INSULATOR 
TRANSITION 

Let us now turn to the Mott metal-insulator transi- 
tion 1 22 3^] from a paramagnetic metal to a paramag- 
netic insulator. This transition is found in various tran- 
sition metal oxides, such as V2O3 doped with Cr [B3|. 



The mechanism driving the Mott-transition is believed 
to be the local Coulomb repulsion U between electrons 
on the same lattice site, although the details of the tran- 
sition should also be influenced by lattice degrees of free- 
dom. The simplest model to investigate the correlation 
driven metal-insulator transition is the one-band Hub- 
bard model M M 



H 



where ^ 



E 



(10) 



c'i a {cia) denote creation (annihilation) operators 
for a fermion on site i and the i,j are the hopping matrix 
elements between site i and j (j37|. Despite its simple 
structure, the solution of this model turns out to be an 
extremely difficult many-body problem. The situation is 
particularly complicated near the metal-insulator transi- 
tion where U and the bandwidth are roughly of the same 
order such that pcrturbative schemes (in U or t) are not 
applicable. 

The existence of a metal-insulator transition in the 
paramagnetic phase (3^] of the half-filled Hubbard model 
has been known since the early work of Hubbard [j34| . 
The details of the transition, however, remained unclear, 
except in the particular case of dimension d = 1, where 
the transition occurs at U = + [^2 39 . Even in the op- 
posite limit of infinite dimensions, where a numerically 
exact solution of the Hubbard model is in principle pos- 
sible, a general consensus concerning the details of the 
transition scenario has not been reached so far. 

Neglecting the transition to an antiferromagnetic 
phase or suppressing it by frustration two coexisting 
solutions are found in DMFT at very low temperatures, 
one insulating and one metallic |4C|j . The transition is 
then of first order, even in the absence of a coupling to 
lattice degrees of freedom. The scenario of a first-order 
transition was first proposed in Refs. J4l| , [42"|| , within cal- 
culations based on the iterated perturbation theory and 
exact diagonalization. It was later confirmed by the NRG 
for T — p5| and quantum Monte-Carlo calculations for 
T > [Ei |l3 . Criticism of this scenario can be found in 
Refs. [||]44p|4^]. 



The results from the NRG for the T = metal-insulator 
transition can be summarized as follows (for details see 
Jl5[). On approaching the transition from the metallic 
side, a typical three-peak structure appears in the spec- 
tral function, with upper and lower Hubbard bands at 
lo w ±t//2 and a quasiparticle peak at lo = 0. The 
width of the quasiparticle peak vanishes for U — > U C 2, 
leaving behind two well-separated Hubbard peaks (see 
Fig. 2 in Q). Although the NRG is not able to resolve 
a small spectral weight between the Hubbard peaks, the 
results indicate that the gap opens discontinuously (see 
also ||). On decreasing U, the transition from the in- 
sulator to the metal occurs at a lower critical value U c \, 
where the gap vanishes. Concerning the numerical value 
of U C 2 ~ 1A7W (W: bandwidth) excellent agreement 
with the result from the projective self-consistent method 
@,| is found. 

The extension of the NRG to T > will now be used 
to determine the full shape of the hysteresis region non- 
perturbatively. 



IV. RESULTS 
A. Spectral function for T > T c 

Figure ^ shows the spectral function A(u>) for vari- 
ous values of U at T = 0.0276W 7 . This is above the 
temperature of the critical point (which we estimate as 
T c sa 0.02VF), so that there is no real transition but a 
crossover from a metallic-like to an insulating-like so- 
lution. As already found in Refs. 0^3], the crossover 
region is nevertheless very narrow, with a very rapid 
suppression of the quasiparticle peak. This is seen also 
in the NRG results (Fig. |^) when U is increased from 
U = 1.05W to U = 1.20VF. The spectral weight of the 
quasiparticle peak is gradually redistributed and shifted 
to the upper (lower) edge of the lower (upper) Hub- 
bard band. An additional structure within the Hubbard 
bands, as reported in |||3|, ^ s not f° un d and would be 
very difficult to see due to the limited resolution of the 
NRG at higher frequencies. 

The inset of Fig. [| shows the Independence of the value 
of the spectral function at zero frequency A(lu = 0). For 
higher values of U, the spectral density at the Fermi level 
is still finite and vanishes only in the limit U — > 00 (or 
for T -> 0, provided that U > U c2 {T = 0)). 

The [/-dependence of A(uj = 0) is shown in Fig. ^a 
for different temperatures. As discussed in Sec. II, the 
temperatures are chosen as T m — T\ ■ A™ 1 , with T\ = 
0.0168VF and m = 0,1,2,3, (A = 1.64 is used for all 
results shown in this paper; the number of states retained 
after truncation is 600, not counting degeneracies). The 
derivative of A(oj = 0) with respect to U, 



A'{uj = 0) 



_ dA(u = 0, U) 
= dU 



(11) 
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is plotted in Fig. f|b. The {/-value where \A'(lu = 0)| 
reaches its maximum defines a characteristic interac- 
tion strength U* for the crossover from metallic-like to 
insulating- like behavior in the region T>T C ; for the def- 
inition of t/ c i,2 for T <T C , see below. Furthermore, the 
width AU of the crossover region can be defined as the 
width at half-height of the peak in A'(oj = 0). 
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FIG. 3. Spectral function for the half-filled Hubbard model 
for various values of U at T — 0.02761T > T c (in the crossover 
region). The crossover from the metal to the insulator occurs 
via a gradual suppression of the quasiparticle peak at ui = 0. 
The inset shows the (7-dependence of A(uu = 0), in particular 
the rapid decrease for U « 1.1 IT. 
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FIG. 4. (a): The [/-dependence of A{uj = 0) for different 
temperatures; the data for T = 0.0168IT < T c show a very 
small hysteresis, not visible on this scale. The other three sets 
of data are for T > T c . (b): The derivative of A(cu = 0,U) 
with respect to U for the same temperatures as in (a). 



Upon lowering the temperature, the width AU rapidly 
decreases and vanishes at the critical temperature T c , 
since A'(u> — 0) diverges at T c (this feature has already 
been discussed in f48|| ). A precise value for T c cannot be 
given as we are presently not able to vary the temperature 
within NRG continuously. The critical temperature is 
estimated as T c w 0.02VT, as a very small hysteresis is 
still present for T = 0.0168VF (on the scale of Fig. |, U cl 
and U C 2 cannot be distinguished). 

The U* as defined above slowly decreases upon increas- 
ing the temperature. This is not at variance with the 
opposite trend observed in Refs. fl6|,f7p3|] (in Ref. |23|], the 
slope of U* changes sign at T » 0.25VF) and depends 
on the definition of U* . Taking U* as, e.g., the value 
of U where A(u> = 0) has dropped to 1% of its value at 
U = would result in an increase of U* upon increasing 
the temperature. 



B. Breakdown of Fermi Liquid versus 
Metal-Insulator Transition 

We now discuss the question of how to define a useful 
criterion for the metal-insulator transition at finite tem- 
peratures. At zero temperature a suitable criterion is the 
vanishing, with increasing U, of the quasiparticle weight 
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duj 



(12) 



The physical meaning of Z is clear for the paramagnetic 
state at T = 0, where the system is either a Fermi liquid 
(for U < U c ) or an insulator (for U > U c ). The vanish- 
ing of Z therefore marks the metal-insulator transition 
at T = 0, as discussed, e.g, in |]|l5|]. This criterion, 
however, cannot be taken over straightforwardly to fi- 
nite temperatures, since for T > the breakdown of the 
Fermi liquid state and the appearence of the insulating 
state do not coincide. Although thispoint has been noted 
before in the literature (see, e.g., |^9| , p3| ), the vanishing 
of Z has been used as (one) criterion for the occurence 
of the metal- insulator transition in p^| . It should be 
noted that the definition of Z used in the finite temper- 
ature quantum Monte-Carlo calculations of Ref. j2ij is 
different from Eq. ( p"2| ) since 9Re£(w)/<%;| w=0 was ap- 
proximated by ImE(io;o)/wo, with o>o the first Matsubara 
frequency. 

To elucidate this point, it is instructive to discuss the 
behavior of the self-energy in the crossover region from 
the metallic-like to the insulating-like solution. The real 
and imaginary part of E(a>) are shown in Fig. 
the same temperature and U values as in Fig. 
U = 1.05W and U — 1A0W the imaginary part shows 
the characteristic structure of the self-energy for a Fermi 
liquid (with the oj 2 dependence for small frequencies and 
the falling off at higher frequencies which leads to a two- 
peak structure) , but with a rapidly increasing scattering 




G 



rate at to = for increasing U. The two-peak structure 
gradually evolves into a structure with a well-pronounced 
peak at lu — characteristic for an insulating solution 
(a vanishing A(tu = 0) would correspond to a (5-function 
in ImS(w)). Note that for J7 > 1.15W the value of 
ImE(cJ = 0) is much larger than the T 2 -term observed 
for T — ► 0. Hence, the mechanism for the strong scat- 
tering at a; = is noi the quasiparticle-interaction but is 
caused by the bare local Coulomb repulsion, which de- 
stroys the Fermi liquid behavior for U > 1.15W. For 
U = 1.15W there is still a narrow dip in ImS(w) at oj = Q 
corresponding to the remnant of the quasiparticle peak 
seen in Fig. [| 
0.0 



3 



-1.0 



-2.0 




U/W=1.05 

- U/W=1.10 

- U/W=1.15 

- U/W=1.20 




FIG. 5. Imaginary part (a) and real part (b) of the 
self-energy for the same temperature (T = 0.0276IF) and 
[/-values as in Fig. ^. The slope of ReE(u) changes sign at 
the same [/-value for which the peak at u = appears in 
ImE(cj). 



more delicate, see, e.g., M), The result for lmS(icj) is 
shown in Fig. [| for the same parameters as in Figs. || 
and ||. The circles indicate the value of lmE(iu> n ) for the 
Matsubara frequencies 



1), n = 0,l,2,. 



(14) 



As the self-energy is defined on the whole imagi- 

nary frequency axis, not only for the Matsubara frequen- 
cies, one can, for instance, check the trivial condition 
ImE(icj) = ImE(w) for u> — > 0. 



0.0 



w 
E 



-1.0 



-2.0 



-3.0 




U/W=1 .05 
U/W=1.10 _ 4 .o 
U/W=1.15 
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0.00 0.05 0.10 0.15 



0.0 0.2 0.4 0.6 0.8 

co/W 
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FIG. 6. Imaginary part of the self-energy, ImE, on the 
imaginary frequency axis for the same parameters as in Figs. ^ 
and |B| The values of ImE for the Matsubara frequencies are 
indicated by the circles. The inset focuses on the change of 
sign of the slope in ImE(icj) for values of U = 1.1W up to 
U = 1.17W (from top to bottom). 



For U = 1.05W and U = 1.10W, the corresponding 
real part of £(w) shows the typical Fermi liquid behavior 
with a negative slope at u> = 0. Upon further increas- 
ing the U, however, the slope of ReE(w) changes sign 
right at the U- value where the peak at to — appears in 
lxsxEi(u>); this is obvious from the Kramers-Kronig trans- 
formation which connects real and imaginary part. Note 
that the l/w-behavior in ReE(w) for larger frequencies is 
not visible on this scale. 

From the full frequency dependence of on the real 
axis one can easily perform the analytic continuation to 
Tt(z) for any value of z in the upper complex plane: 



E(z) 



1 t^^J) 



(13) 



In particular for z = iui and ui real, Eq. ( |l3| ) gives the 
real and imaginary parts of the self-energy on the imag- 
inary frequency axis (the analytic continuation from the 
imaginary to the real frequency axis, however, is much 



Furthermore, the slope of ImE(iw) for u> — *• is iden- 
tical to the slope of the real part of E(w). As a conse- 
quence, the same change in the slope of the self-energy 
is visible in both Figs. [|b and || The inset of Fig. g 
illustrates this for a smaller frequency range and a nar- 
row mesh of [/-values from U = 1.1W up to U = 1.17W 
(from top to bottom). 

This behavior of the self-energy has drastic conse- 
quences for the notion of a quasiparticle weight Z in the 
crossover regime. We see that the application of Eq. jl^ ) 
to the self-energies as obtained in Figs. ^| and || leads to 
unphysical results for U > 1.1514^. Due to the change 
of sign in 8R ° S M upon increasing U, the Z from 

Eq. (|12j) starts increasing again and even diverges at a 
particular value of U for which the derivative of the self- 
energy is equal to one. For larger values of U, Z becomes 
negative and approaches zero from below for U — > oo. 
Apparently the use of Eq. ( p^ ) does not make sense for 
U > 1.15W [BO! which is due to the fact the the concept of 
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quasiparticles itself breaks down in the crossover regime. 
The quasiparticle weight is therefore not an appropriate 
measure for the metal-insulator transition in the whole 
parameter space. Note also that in the crossover region, 
the weight of the remnant of the quasiparticle peak is not 
associated to Z. 

Whereas there is no unique criterion for a characteristic 
value U* for T > T c , critical values for U can nevertheless 
be defined for T < T c via the value of U at which A(u> = 
0) changes discontinuously. 



C. Spectral function for T <T C 

Figure ^ shows the spectral function A(ui) in the hys- 
teresis region for T = 0.0103VF, both for increasing U 
(Fig. 0a) and decreasing U (Fig. 0b). The results are 
shown for a very fine mesh of [/-values close to U C 2 ~ 
1.21W and U cl » l.UW. 
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FIG. 7. Spectral function for T = 0.0103W; (a): increasing 
U (b): decreasing U; the transitions at U C 2 ~ 1-21W and 
Uci ~ 1.14 W are characterized by a significant redistribution 
of spectral weight and a jump in A(u> = 0) (see also Fig. 0) . 

In both cases, the transition is of first order, i.e., as- 
sociated with a discontinuous redistribution of spectral 
weight. The hysteresis effect is further illustrated in the 
[/-dependence of A(u = 0) for T = 0.01031F < T c (Fig.|). 

Whereas the critical values U c i and [Z c2 can be easily 
defined by the jump in A(ui = 0), the calculation of the 
actual thermodynamic transition requires the knowledge 
of the free energy F of both metallic and insulating so- 
lutions. The determination of F goes beyond the scope 
of this paper. There is no way of directly calculating F 
within the NRG approach, so one has to determine the 
free energy via integrating over a path from a particu- 
lar point in the (U, T)-plane for which the free energy 



is known, up to the actual values of U and T. However, 
the knowledge of the U c (T) for the actual thermodynamic 
transition will not alter the fact that the transition is of 
first order. 



1.5 



1.0 



< 0.5 



0.0 
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A V 



0.4 



0.8 



1.2 



1.6 



U/W 



FIG. 8. [/-dependence of A(lo = 0) for T = 0.0103W; solid 
line: increasing U, dashed line: decreasing U. 



D. Phase diagram 

Let us finally discuss the phase diagram for the Mott 
metal-insulator transition in the very low temperature 
region. 

In Fig. ||, the dashed lines for T > T c indicate the posi- 
tion and width of the crossover region as calculated from 
the NRG data of Fig. ||. The open circles and squares 
are the NRG-results for U C \{T) and U C 2{T), respectively. 
As the NRG calculations cannot, so far, be performed for 
arbitrary values of T, we cannot give a precise value for 
the critical point. The U C 2{T) nicely extrapolates to the 
previously obtained value for T = jl5j ; the same is true 
for U cl {T). Note that the value for U cl (T = 0) = 1.195W 
plotted here is slighly reduced as compared to the orig- 
inally published value [/ cl (T = 0) = 1.25W ||. This is 
due to the different value for A, the number of states and 
the broadening used here. 

Figure ^ also contains recent quantum Monte-Carlo 
results of Joo and Oudovenko [Q (filled symbols), as 
well as the result from the iterated perturbation theory 
H which tends to overestimate both U c i(T) and U C 2(T). 
The phase boundaries obtained from the NRG are below 
the values obtained from the quantum Monte-Carlo re- 
sults. Concerning the NRG values, it is well-known that 
due to the logarithmic discretization the NRG tends to 
underestimate the effective hybridization [5])] (hence un- 
derestimating the value of U necessary to overcome the 
kinetic energy). This effect has, e.g., been studied in the 
context of the quantum phase-transition from the local 
moment to the strong coupling phase in the soft-gap An- 
derson model [fL9|. For the transition in this model, the 
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value of U c for A = 2.0 is about 5% below the extrapo- 
lated value for A — ► 1; more importantly, the U C (A) is a 
perfectly straight line from A=1.4 to A = 3.0. 
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FIG. 9. Results for the phase diagram of the Mott transi- 
tion obtained from different methods: NRG (open symbols), 
quantum Monte-Carlo (QMC, filled symbols) and iterated 
perturbation theory (IPT, solid lines). The dashed lines for 
T > T c indicate the position and width of the crossover region 
as calculated from the data of Fig. ^. 

A similar A — > 1 extrapolation is difficult to perform 
for the metal-insulator transition studied here, since al- 
ready a very large number of DMFT iterations is neces- 
sary to determine a single value of U C (A). Calculations of 
U c i and U C 2 for one value of T with A=1.64 and A = 2.0 
at least show the expected trend, i.e., a slight increase of 
the U c 's with decreasing A. 

Taking into account the unavoidable numerical errors 
in both procedures, the agreement between NRG and 
quantum Monte-Carlo results for the phase-boundary is 
seen to be very good; the agreement can even be further 
improved Q. 

V. SUMMARY 

In this paper we presented results from the numeri- 
cal renormalization group method (NRG) for the finite 
temperature Mott transition in the Hubbard model on 
a Bcthe lattice within dynamical mean held theory. For 
the crossover region T > T c , the quasiparticle peak in 
the spectral function gradually vanishes upon increasing 
U and the imaginary part of the self-energy develops a 
sharp peak at lj — 0. Associated with this is a change 
of sign of ReS(w) at u> = 0. As a consequence, the be- 
haviour of the quasiparticle weight can no longer be used 
as a criterion for the transition at finite temperature. 



For T < T c w 0.021^, we find two coexisting solutions 
in the range U cX (T) < U < U c2 {T). The values for the 
critical U can be determined for arbitrarily small temper- 
atures, in contrast to the quantum Monte-Carlo method 
which is so far restricted to T > W/150. The critical 
U c i(T) and U C 2(T) are characterized by a redistribution 
of finite spectral weight in the spectral function. 

We therefore obtain a consistent picture for the Mott 
metal-insulator transition from a paramagnetic metal to 
a paramagnetic insulator in the whole parameter regime. 
The results are in very good agreement with those from 
other non-perturbative methods (the quantum Monte- 
Carlo method and the projective self-consistent method) 
in their respective ranges of applicability. 

There are still several questions left for further inves- 
tigations. A continuous variation of the temperature 
within the NRG requires a better understanding of the 
A-dependence of the results. The NRG also allows the 
calculation of a variety of dynamic and transport prop- 
erties in the whole parameter regime, such as dynamic 
susceptibility and optical conductivity. A generalization 
of the NRG method to antiferromagnetic phases and the 
Hubbard model away from half-filling is in progress. 
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